1. Compositacion de Sondajes
Una variable regionalizada \(Z(\mathbf{u})\) puede definirse, no sólo en cada punto del espacio \(\mathbf{u}\), sino que también en una superficie (2D) o en un volumen (3D), \(Z(\mathbf{V})\). La superficie o el volumen sobre el cual se considera la variable regionalizada se denomina soporte. En general, el soporte de las mediciones es muy pequeño (asimilado a un «punto»), mientras que el que interesa en la práctica puede ser más voluminoso (por ejemplo, las unidades selectivas de explotación en evaluación minera o las unidades de remediación en contaminación de suelo). Esta noción es esencial debido a la dependencia que existe entre el soporte y la distribución estadística de los valores, conocida como efecto soporte: los soportes voluminosos presentan una menor cantidad de valores extremos y una mayor cantidad de valores intermedios que los soportes puntuales. Así, la distribución de los valores (en especial, su varianza) depende del soporte sobre el cual está definida la variable regionalizada.
1.1. Instalar PyGeostat
Instalamos pygeostat, un paquete de Python para modelado geoestadístico. Pygeostat está orientado a la preparación de datos espaciales, la automatización de flujos de trabajo geoestadísticos y el modelado utilizando herramientas desarrolladas en el Centre for Computational Geostatistics (CCG)
[ ]:
## La instalación pide re-iniciar el kernel, asi que no asustarse por esto.
!pip install pygeostat
Desafortunadamente, la librería no ha sido actualizada desde el 2015, por lo que debemos descargar y reemplazar en nuestro computador el archivo desurvey.py. Si estas usando Jupyter, deberas encontra este archivo por tu cuenta en tu computador para reemplazarlo a mano (copiar y pegar el nuevo archivo descargado en tu carpeta).
Si estas usando colab, basta con correl el siguiente código:
[ ]:
!wget -q https://raw.githubusercontent.com/arqlm/arqlm.github.io/main/_static/desurvey.py
!mv desurvey.py /usr/local/lib/python3.12/dist-packages/pygeostat/datautils/desurvey.py
1.2. Importar librerias
[1]:
import pygeostat as gs
import numpy as np
import pandas as pd
[2]:
# Leer collares
collar = gs.DataFile(flname='./dataEvYac_collar.csv',readfl=True,x='SURV_X',
y='SURV_Y',z='SURV_Z',dh='HOLEID')
print('Collars:\n',collar.data.head())
# Leer Surveys
survey = gs.DataFile(flname='./dataEvYac_survey.csv',readfl=True,dh='HOLEID')
print('\nSurveys:\n',survey.data.head())
# Leer Assays: puedes usar tu data una vez que hallas encontrado las ug's
rawdata = gs.DataFile(flname='./Data_sin_compositar_con_UGs.csv',readfl=True,dh='dhid',
ifrom='from',ito='to')
print('\nRaw data file:\n',rawdata.data.head())
rawdata.data = rawdata.data[rawdata.data['cu_pct'] >= 0] # Filtrar valores negativos
rawdata.data
Collars:
HOLEID SURV_X SURV_Y SURV_Z DEPTH
0 CCDDH-001 472560.798 6925645.753 4204.009 800.10
1 CCDDH-002A 472479.288 6925449.980 4289.060 1010.30
2 CCDDH-003 471898.517 6925719.817 4287.854 860.22
3 CCDDH-004 472109.838 6925708.225 4278.548 1100.00
4 CCDDH-005 472233.539 6925602.721 4322.216 918.00
Surveys:
HOLEID Along Azimuth Inclination
0 CCDDH-001 0.0 209.000 -60.326
1 CCDDH-001 10.0 208.094 -60.150
2 CCDDH-001 20.0 208.210 -60.143
3 CCDDH-001 30.0 208.512 -60.252
4 CCDDH-001 40.0 208.276 -60.215
Raw data file:
dhid length from to Este Norte Elevación au_ppm \
0 98CCD089 12.2 0.0 12.2 472186.686 6925804.447 4220.763 -99.00
1 98CCD089 1.8 12.2 14.0 472187.202 6925805.493 4213.861 0.30
2 98CCD089 2.0 14.0 16.0 472187.343 6925805.770 4211.986 0.47
3 98CCD089 2.0 16.0 18.0 472187.493 6925806.060 4210.013 0.31
4 98CCD089 2.0 18.0 20.0 472187.642 6925806.351 4208.040 0.29
ag_ppm cu_pct aucn_ppm cucn_ppm Zmin Alte Lito UG
0 -99.0 NaN -99.0 -99.0 OXI ILL_CLO IBX_MM 2.0
1 2.3 0.011 -99.0 -99.0 OXI ILL_CLO IBX_MM 2.0
2 16.2 0.032 -99.0 -99.0 OXI ILL_CLO IBX_MM 2.0
3 2.3 0.018 -99.0 -99.0 OXI ILL_CLO IBX_MM 2.0
4 2.1 0.010 -99.0 -99.0 OXI ILL_CLO IBX_MM 2.0
[2]:
| dhid | length | from | to | Este | Norte | Elevación | au_ppm | ag_ppm | cu_pct | aucn_ppm | cucn_ppm | Zmin | Alte | Lito | UG | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 98CCD089 | 1.80 | 12.2 | 14.00 | 472187.202 | 6925805.493 | 4213.861 | 0.30 | 2.3 | 0.011 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 2.0 |
| 2 | 98CCD089 | 2.00 | 14.0 | 16.00 | 472187.343 | 6925805.770 | 4211.986 | 0.47 | 16.2 | 0.032 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 2.0 |
| 3 | 98CCD089 | 2.00 | 16.0 | 18.00 | 472187.493 | 6925806.060 | 4210.013 | 0.31 | 2.3 | 0.018 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 2.0 |
| 4 | 98CCD089 | 2.00 | 18.0 | 20.00 | 472187.642 | 6925806.351 | 4208.040 | 0.29 | 2.1 | 0.010 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 2.0 |
| 5 | 98CCD089 | 2.00 | 20.0 | 22.00 | 472187.785 | 6925806.636 | 4206.066 | 0.40 | 2.1 | 0.010 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 2.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 77836 | ZVD002 | 2.00 | 684.0 | 686.00 | 471293.198 | 6925823.101 | 3761.674 | 0.03 | 1.0 | 0.002 | -99.0 | -99.0 | BACK | PROP | ECS | 2.0 |
| 77837 | ZVD002 | 2.00 | 686.0 | 688.00 | 471293.540 | 6925824.040 | 3759.942 | 0.01 | 0.7 | 0.001 | -99.0 | -99.0 | BACK | PROP | ECS | 2.0 |
| 77838 | ZVD002 | 2.00 | 688.0 | 690.00 | 471293.882 | 6925824.980 | 3758.209 | 0.01 | 0.4 | 0.001 | -99.0 | -99.0 | BACK | PROP | ECS | 2.0 |
| 77839 | ZVD002 | 2.00 | 690.0 | 692.00 | 471294.224 | 6925825.920 | 3756.477 | 0.06 | 1.3 | 0.004 | -99.0 | -99.0 | BACK | PROP | ECS | 2.0 |
| 77840 | ZVD002 | 2.48 | 692.0 | 694.48 | 471294.607 | 6925826.972 | 3754.538 | 0.06 | 0.5 | 0.003 | -99.0 | -99.0 | BACK | PROP | ECS | 2.0 |
76907 rows × 16 columns
Ahora configuramos el largo del compósito final a obtener: en este caso, 6 metros:
[3]:
%%capture
## muteamos la salida para que no sea tan larga (puede borrar %%capture si quieres ver todo)
comps = gs.desurvey.set_comps(rawdata,6.0) ## Segundo parámetro controla el largo del compósito final
upscaled = gs.desurvey.get_comps(comps,rawdata,{'cu_pct':'continuous','UG':'categorical'}) ## El tercer parámetro indica el tipo de variable (continua o categórica)
rawdata.origdata = rawdata.data.copy() # Guardar una copia de los datos originales por si acaso...
rawdata.data = upscaled # Hacer que los datos ampliados sean el conjunto de datos principal
Ya podemos ver las estadisticas de la data en su nuevo volumen:
[4]:
upscaled_pos = upscaled[upscaled['cu_pct'] >= 0] # Filtrar valores negativos
upscaled_pos.describe()
[4]:
| cu_pct | UG | from | to | |
|---|---|---|---|---|
| count | 25414.000000 | 25414.000000 | 25414.000000 | 25414.000000 |
| mean | 0.145067 | 2.944519 | 358.251330 | 364.251330 |
| std | 0.183558 | 0.999070 | 296.482491 | 296.482491 |
| min | 0.001000 | 1.000000 | 0.000000 | 6.000000 |
| 25% | 0.022667 | 2.000000 | 108.000000 | 114.000000 |
| 50% | 0.086454 | 3.000000 | 272.000000 | 278.000000 |
| 75% | 0.219667 | 4.000000 | 564.000000 | 570.000000 |
| max | 6.723333 | 4.000000 | 1470.000000 | 1476.000000 |
Solo falta agregar la posicion o centroide a cada muestra compositada:
[5]:
drillholes = gs.desurvey.set_desurvey(collar,survey,'Along','Azimuth','Inclination')
WARNING: drillhole routines are untested with recent Pygeostat changes!
Proceed with caution!
[6]:
drillhole_items = list(drillholes.items())
[7]:
%%capture
## muteamos la salida para que no sea tan larga (puede borrar %%capture si quieres ver todo)
gs.desurvey.get_desurvey(rawdata,drillholes,x=collar.x,y=collar.y,z=collar.z)
print(rawdata.data.head())
Filtramos algunos algunos tramos de sondaje sin collares:
[8]:
df_filtered = rawdata.data[rawdata.data['SURV_X'] > -100] # Filtrar valores no deseados
df_filtered
[8]:
| cu_pct | UG | dhid | from | to | SURV_X | SURV_Y | SURV_Z | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0.020300 | 2.0 | 98CCD089 | 12.2 | 18.2 | 472187.343570 | 6.925806e+06 | 4211.772172 |
| 1 | 0.009500 | 2.0 | 98CCD089 | 18.2 | 24.2 | 472187.740701 | 6.925807e+06 | 4205.844155 |
| 2 | 0.006767 | 2.0 | 98CCD089 | 24.2 | 30.2 | 472188.089349 | 6.925807e+06 | 4199.915106 |
| 3 | 0.010433 | 2.0 | 98CCD089 | 30.2 | 36.2 | 472188.470072 | 6.925808e+06 | 4193.995578 |
| 4 | 0.011133 | 2.0 | 98CCD089 | 36.2 | 42.2 | 472188.878860 | 6.925809e+06 | 4188.084383 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 22343 | -999.000000 | -999.0 | CCDDH-007 | 558.0 | 564.0 | 471959.484334 | 6.925444e+06 | 3935.475961 |
| 22344 | -999.000000 | -999.0 | CCDDH-007 | 564.0 | 570.0 | 471961.631627 | 6.925447e+06 | 3930.686670 |
| 22345 | -999.000000 | -999.0 | CCDDH-007 | 570.0 | 576.0 | 471963.789016 | 6.925450e+06 | 3925.900536 |
| 22346 | -999.000000 | -999.0 | CCDDH-007 | 576.0 | 582.0 | 471965.956499 | 6.925453e+06 | 3921.117560 |
| 22347 | -999.000000 | -999.0 | CCDDH-007 | 582.0 | 588.0 | 471968.165175 | 6.925455e+06 | 3916.340916 |
16927 rows × 8 columns
Y finalmente, filtramos algunos algunos tramos de sondaje sin ensayos:
[9]:
upscaled_pos = df_filtered[df_filtered['cu_pct'] >= 0] # Filtrar valores negativos
upscaled_pos.describe()
[9]:
| cu_pct | UG | from | to | SURV_X | SURV_Y | SURV_Z | |
|---|---|---|---|---|---|---|---|
| count | 16882.000000 | 16882.000000 | 16882.000000 | 16882.000000 | 16882.000000 | 1.688200e+04 | 16882.000000 |
| mean | 0.165559 | 3.143822 | 391.293660 | 397.293660 | 472141.581163 | 6.925472e+06 | 3966.341566 |
| std | 0.196619 | 0.926312 | 284.212204 | 284.212204 | 302.008871 | 2.401475e+02 | 258.580610 |
| min | 0.001000 | 2.000000 | 0.000000 | 6.000000 | 471116.243127 | 6.924722e+06 | 2950.021900 |
| 25% | 0.032333 | 2.000000 | 153.190000 | 159.190000 | 471927.824777 | 6.925302e+06 | 3787.725996 |
| 50% | 0.116437 | 4.000000 | 330.310000 | 336.310000 | 472169.814450 | 6.925478e+06 | 4011.448263 |
| 75% | 0.242961 | 4.000000 | 597.475000 | 603.475000 | 472356.513996 | 6.925649e+06 | 4171.006788 |
| max | 6.723333 | 4.000000 | 1401.050000 | 1407.050000 | 473023.885835 | 6.926165e+06 | 4445.782649 |
Es posible ahora visualizar la información compositada en el espacio:
[10]:
import plotly.io as pio
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py
[11]:
df_filtered["UG"] = df_filtered["UG"].astype(str)
fig = px.scatter_3d(df_filtered, x='SURV_X', y='SURV_Y', z='SURV_Z', color='UG',category_orders={'UG': sorted(df_filtered.groupby('UG').groups.keys())})
fig.update_traces(marker=dict(size=1.0))
# Mejor estilo para leyenda
fig.update_layout(
title_text='Unidades Geológicas para Cu',
legend_title_text='UGs',
legend=dict(
itemsizing='constant',
font=dict(size=14),
bgcolor='rgba(255,255,255,0.8)',
bordercolor='black',
borderwidth=1
)
)
fig.show()
C:\Users\Alvaro\AppData\Local\Temp\ipykernel_17776\2927694018.py:1: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
Data type cannot be displayed: application/vnd.plotly.v1+json
[12]:
continuous = 'cu_pct'
import plotly.io as pio
# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py
fig = go.Figure(data=[go.Scatter3d(
x=df_filtered['SURV_X'],
y=df_filtered['SURV_Y'],
z=df_filtered['SURV_Z'],
marker=dict(color=df_filtered[continuous],
colorscale=px.colors.sequential.Rainbow[1:],
cmin=0.0,
cmax=df_filtered[continuous].quantile(0.95),
size=1.0, # Esta opcion controla el tamaño de los datos
colorbar=dict(
title='Cu [%]',
thickness=20,
# Add a border to the colorbar
outlinecolor='black', # Color of the border
outlinewidth=2, # Width of the border
bordercolor='white', # Background border color (if applicable)
borderwidth=1 # Background border width
)
),
mode='markers',
opacity=1
)])
# otras opciones para controlar el color de los elementos de la figura
fig.update_layout(scene = dict(
xaxis = dict(
title='Este',
backgroundcolor="white",
gridcolor="gray",
showbackground=True,
zerolinecolor="white",),
yaxis = dict(
title='Norte',
backgroundcolor="white",
gridcolor="gray",
showbackground=True,
zerolinecolor="white"),
zaxis = dict(
title='Elevación',
backgroundcolor="white",
gridcolor="gray",
showbackground=True,
zerolinecolor="white",),
),
font_family="Times New Roman",
font_color="black", hovermode=False
)
fig.show()
[13]:
rawdata.data
[13]:
| cu_pct | UG | dhid | from | to | SURV_X | SURV_Y | SURV_Z | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0.020300 | 2.0 | 98CCD089 | 12.20 | 18.20 | 472187.343570 | 6.925806e+06 | 4211.772172 |
| 1 | 0.009500 | 2.0 | 98CCD089 | 18.20 | 24.20 | 472187.740701 | 6.925807e+06 | 4205.844155 |
| 2 | 0.006767 | 2.0 | 98CCD089 | 24.20 | 30.20 | 472188.089349 | 6.925807e+06 | 4199.915106 |
| 3 | 0.010433 | 2.0 | 98CCD089 | 30.20 | 36.20 | 472188.470072 | 6.925808e+06 | 4193.995578 |
| 4 | 0.011133 | 2.0 | 98CCD089 | 36.20 | 42.20 | 472188.878860 | 6.925809e+06 | 4188.084383 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 25467 | 0.006253 | 2.0 | ZVD002 | 665.28 | 671.28 | -999.000000 | -9.990000e+02 | -999.000000 |
| 25468 | 0.017613 | 2.0 | ZVD002 | 671.28 | 677.28 | -999.000000 | -9.990000e+02 | -999.000000 |
| 25469 | 0.001000 | 2.0 | ZVD002 | 677.28 | 683.28 | -999.000000 | -9.990000e+02 | -999.000000 |
| 25470 | 0.001333 | 2.0 | ZVD002 | 683.28 | 689.28 | -999.000000 | -9.990000e+02 | -999.000000 |
| 25471 | 0.003108 | 2.0 | ZVD002 | 689.28 | 695.28 | -999.000000 | -9.990000e+02 | -999.000000 |
25472 rows × 8 columns
Ya sabes como exportar la información en un archivo CSV y como hacer un histograma, para verificar los efectos del nuevo soporte en la distribución.
[14]:
upscaled_pos.to_csv('Data_compositada_con_UGs.csv',index=False)